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ABSTRACT 

A  major  computational  code  called  CERV  was 
developed  to  determine  complex  equilibrium  compo¬ 
sitions  of  a  nonideal  mixture  of  numerous  imperfect 
gases  and  compressible  liquid  and  solid  species  with 
phase  transitions  for  closed-vessel  applications.  This 
code  minimizes  Gibbs  energy  using  reaction  variables, 
in  contrast  to  other  major  codes  like  BRL-Blake,  BRC- 
Bagheera  and  NASA-CEA  that  use  composition  vari¬ 
ables  such  as  mole  numbers.  The  CERV  code  has 
significant  advantages  in  handling  compressible  con¬ 
densed  species  with  phase  transitions  and  comput¬ 
ing  nonideal  mixture  compositions  accurately  and  effi¬ 
ciently.  Computations  are  done  robustly,  without  con¬ 
vergence  failures  from  matrix  inversions  or  iterative 
procedures,  for  problems  with  reaction  products  con¬ 
sisting  of  hundreds  of  gaseous,  liquid  and  solid  species. 

1.  INTRODUCTION 

A  thorough  knowledge  of  chemical  equilibrium 
mixture  compositions  and  thermodynamic  properties 
of  the  reaction  products  of  energetic  materials  is  im¬ 
portant  not  only  in  the  design,  performance  analy¬ 
sis  and  utilization  of  engineering  equipment  such  as 
compressors,  turbines,  heat  exchangers,  combustion 
engines,  shock  tubes,  rockets  and  projectile  launch¬ 
ers,  but  also  in  the  formulation  and  characteriza¬ 
tion  of  energetic  materials  such  as  fuel-air  mixtures, 
propellants,  explosives,  primers,  fuses  and  igniters. 
The  equilibrium-mixture  composition  (i.e.,  types  of 
gaseous,  liquid  and  solid  species  in  the  reaction  prod¬ 
ucts  and  their  mole  numbers)  and  its  thermodynamic 
properties  (e.g.,  pressure,  volume,  temperature,  in¬ 
ternal  energy,  Gibbs  energy)  are  normally  not  avail¬ 
able  from  measurements  or  experiments,  so  they  are 
predicted  instead  by  using  computational  methods 
based  on  appropriate  models  for  the  chemical  activ¬ 
ity,  type  of  mixture  (ideal  or  nonideal)  and  pressure- 
volume-temperature  properties  (equation  of  state)  of 
the  gaseous,  liquid  and  solid  species.  The  benefits  of 
the  computational  methods  in  providing  accurate  pre¬ 
dictions  for  the  applications  mentioned  previously  de¬ 
pend  considerably  on  the  physical  validity  and  com¬ 
pleteness  of  the  underlying  models. 


The  challenges  of  providing  a  thorough  knowledge 
of  chemical  equilibrium  mixture  compositions,  thermo¬ 
dynamic  properties  and  models  are  being  met  by  ad¬ 
vancements  in  chemical  equilibrium  theory,  modelling 
and  solution  methods  and  by  developments  in  com¬ 
puter  algorithms.  Although  ideal  and  nonideal  mix¬ 
ture  models  are  well  developed,  advancements  are  re¬ 
quired  for  more  accurate  equations  of  state  and  ad¬ 
ditional  thermochemical  data  for  imperfect  gas  be¬ 
haviour  and  for  the  compressibility  of  liquids  and 
solids.  This  paper  illustrates  the  capabilities  and  ad¬ 
vantages  of  a  new  computational  code  (CERV)  in  solv¬ 
ing  combustion  problems  involving  energetic  materi¬ 
als  by  using  present-day  models  for  ideal  and  nonideal 
mixtures,  for  equations  of  state  for  perfect  and  imper¬ 
fect  gases,  and  for  liquid-solid  and  solid-solid  phase 
transitions. 

2.  MAJOR  COMPUTATIONAL  CODES 

The  NASA-CEA  code  was  developed  at  the  NASA 
Lewis  Research  Center  in  the  United  States  (US)  over 
the  past  50  years  to  solve  combustion-product  compo¬ 
sitions  and  thermodynamic  properties  associated  with 
energetic  materials,  rocket  engines  and  shock  waves 
and  detonations  in  gases  for  the  US  aerospace  and  mil¬ 
itary  programs.  This  code  is  based  on  an  ideal  mixture 
model,  an  equation  of  state  (EOS)  for  perfect  gases, 
and  condensed  species  that  have  no  volume  in  the  re¬ 
action  products.  The  Blake  code  was  developed  in 
the  1970s  and  1980s  at  the  Ballistic  Research  Labora¬ 
tory  of  the  US  Army  for  military  applications  involv¬ 
ing  propellants,  explosives,  primers,  fuses  and  igniters. 
Because  the  pressures  and  densities  of  the  reactant- 
product  mixtures  for  such  applications  are  generally 
high,  the  Blake  code  is  based  on  a  nonideal  mixture 
model,  an  EOS  for  imperfect  gases  (e.g.,  virial  EOS), 
and  compressible  condensed  species  that  have  finite 
volumes.  The  Bagheera  code  was  developed  in  the 
1980s  and  1990s  in  France  for  similar  purposes,  and  it 
is  based  on  an  ideal  mixture  model,  a  modified  virial 
EOS  and  includes  the  volumes  of  condensed  species. 
By  the  early  1990s  the  Bagheera  code  supplanted  the 
Blake  code  as  the  standard  for  military  applications 
in  NATO  countries.  See  Wong,  Gottlieb  and  Fussier 
(2003)  for  a  more  detailed  review  of  these  codes. 
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The  CERV  code  was  developed  in  the  1990s  at  the 
University  of  Toronto  in  Canada,  starting  with  the 
doctoral  studies  of  F.C.H.  Wong  (2001),  for  solving 
fairly  general  chemical  equilibrium  problems,  but  with 
emphasis  on  energetic  material  combustion  in  closed 
vessels  for  military  applications.  This  code  has  the 
capability  of  using  either  ideal  or  nonideal  mixture 
models,  different  EOS  for  perfect  and  imperfect  gases 
and  either  incompressible  or  compressible  condensed 
species.  This  code  is  based  on  the  earlier  work  of  Vil- 
lars  (1959,  1960),  Cruise  (1964),  Smith  and  Missen 
(1967,  1968),  who  solved  problems  for  a  specified  mix¬ 
ture  temperature  and  pressure  and  using  the  simplest 
ideal  mixture  of  perfect  gases  and  condensed  species 
without  volume.  Smith  and  Missen  used  the  name 
VCS  method;  we  use  the  acronym  CERV  stemming 
from  chemical  equilibrium  using  reaction  variables. 

The  essential  capabilities  of  the  NASA-CEA, 
Blake,  Bagheera  and  CERV  codes  are  summarized  in 
Table  1.  Most  of  these  codes  can  solve  the  two  basic 
problems  of  a  specified  equilibrium-mixture  tempera¬ 
ture  and  pressure  {T-P  problem)  and  a  specified  en¬ 
ergy  and  volume  {E-V  or  closed- vessel  problem).  If 
individual  volumes  of  condensed  species  are  ignored, 
then  the  modelling  of  condensed  species  as  incom¬ 
pressible  or  compressible  is  irrelevant  or  not  applicable 
(NA) .  Some  codes  are  used  with  well-known  equations 
of  state  stemming  from  Redlich,  Kwong  and  Soave 
(RKS);  Benedict,  Webb,  Rubin,  Starling  and  Han 
(BWRSH);  Becker,  Kistiakowski  and  Wilson  (BKW); 


Table  1.  Capabilities  of  major  computational  codes. 


feature  \  code 

CEA 

Blake 

Bagheera 

CERV 

T-P  problem 

yes 

yes 

no 

yes 

E-V  problem 

yes 

yes 

yes 

yes 

ideal  mixture, 
perfect  gases 

yes^ 

yes^ 

yes^ 

yes^ 

ideal  mixture, 
imperfect  gases 

no 

no 

yes^ 

yes^’^ 

nonideal  mixture, 
imperfect  gases 

no 

4 

yes 

no 

yes^ 

one  phase  for  all 
gaseous  species 

yes 

yes 

yes 

yes 

one  phase  for  each 
condensed  species 

yes 

yes 

yes 

yes 

volumes  of 
condensed  species 

no 

yes 

yes 

yes 

incompressible 
condensed  species 

NA 

yes 

yes 

yes 

compressible 
condensed  species 

NA 

yes 

no 

yes 

phase  transitions 

yes 

yes 

yes 

yes 

different  densities 
for  two  species  in 
phase  transition 

NA 

yes 

no 

yes 

^EOS  for  perfect  gases  ^^irial,  RKS,  BWRSH  EOS 
^modified  virial  EOS  ^virial,  BWK,  PWHK  EOS 


and  Powell,  Wilmot,  Haar  and  Klein  (PWHK).  See 
Wong,  Gottlieb  and  Lussier  (2003)  for  details  and  ref¬ 
erences  on  these  EOS  and  the  virial  EOS. 

3.  BNR  AND  CERV  SOLUTION  METHODS 

The  NASA-CEA,  Blake  and  Bagheera  codes  de¬ 
termine  equilibrium-mixture  compositions  for  T-P  and 
E-V  problems  by  minimizing  the  Gibbs  or  Helmholtz 
energy  of  a  mixture  using  composition  variables  (mole 
numbers).  The  Gibbs  energy  G(n),  where  n  denotes 
all  gaseous  and  condensed  species  mole  numbers  (num¬ 
bering  Ns),  is  minimized  subject  to  mass  conservation 
constraints  on  atomic  elements  that  are  included  ex¬ 
plicitly  by  means  of  a  set  of  Lagrange  constraint  mul¬ 
tipliers  A,  one  for  each  atomic  element  (numbering 
Nf.)  in  the  reactants.  This  results  in  the  Lagrangian 

L(n,  A)  =  G{n)  -  YhZi  “  Y/j=i  >  where 

are  the  entries  of  an  atomic  element  abundance  vec¬ 
tor  a,  and  Fij  are  the  entries  of  a  molecular  formula 
matrix  F.  The  minimization  of  L{n,  A)  yields  the  min¬ 
imum  of  G{n)  and  the  equilibrium  composition  n. 

Several  solution  methods  have  been  developed  for 
minimizing  L{n,\).  The  most  renown  are  the  Brink- 
ley,  NASA  and  RAND  Corporation  methods,  referred 
to  collectively  as  the  BNR  methods,  which  stem  re¬ 
spectively  from  the  work  of  Brinkley  (1947);  Huff, 
Gordon  and  Morrell  (1951);  and  White,  Johnson  and 
Dantzig  (1958).  In  BNR  methods,  an  appropriate  set 
of  gaseous,  liquid  and  solid  species  is  selected  as  reac¬ 
tion  products,  their  mole  numbers  are  estimated,  and 
an  iterative  procedure  (second-order  Newton-Raphson 
scheme)  is  used  to  descend  step  by  step  to  the  mini¬ 
mum  of  L{n,  A).  At  each  step  along  the  descent  path 
a  set  of  linear  equations  of  the  form  AAx  —  b  must  be 
solved  to  yield  Aa;  =  A~^b,  where  Aa;  denotes  correc¬ 
tions  to  both  n  and  A  from  the  previous  step.  Matrix 
A  must  be  invertible  (nonsingular)  or  a  convergence 
failure  will  occur  before  Ax  can  be  evaluated. 

The  minimization  problem  for  BNR  methods  must 
be  posed  carefully  to  avoid  theoretical  singularities  in 
matrix  A.  For  example,  the  initial  selection  of  gaseous 
and  condensed  species  must  not  include  two  condensed 
species  with  the  same  molecular  formula  such  as  liquid 
water  and  ice,  and  the  phase  rule  must  be  satisfied  by 
restricting  the  number  of  condensed  species  in  the  re¬ 
action  products  (normally  to  one  less  than  the  number 
of  atomic  elements  Ng).  Even  when  matrix  A  is  the¬ 
oretically  invertible  (nonsingular),  it  can  become  ill- 
conditioned  or  numerically  singular  during  computa¬ 
tions  from  number  round-off  errors,  and  thereby  result 
in  a  convergence  failure.  This  situation  worsens  dra¬ 
matically  as  the  mass  content  of  the  gaseous  species  di¬ 
minishes  significantly  relative  to  that  of  the  condensed 
species,  a  frequent  occurrence  in  solving  T-P  and  E-V 
problems  involving  primers,  fuses  and  igniters. 

The  most  advanced  and  popular  BNR  methods 
described  by  Gordon  and  McBride  (1994)  contain  a 
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multiple-loop  iterative  scheme,  which  is  embedded  in 
the  NASA-CEA  code  and  emulated  in  the  Bagheera 
code.  To  avoid  theoretical  matrix  singularities  and 
phase-rule  violations,  the  initial  reaction  products  for 
T-P  and  E-V  problems  are  selected  as  gaseous  species 
only  and  their  mole  numbers  are  estimated  (e.g.,  each 
is  one  mole  divided  by  the  number  of  gaseous  species). 
L(n,  A)  and  thereby  G(n)  are  minimized  iteratively  to 
determine  the  corresponding  equilibrium  composition 
n,  which  is  the  first  iterative  loop.  One  condensed 
species  is  possibly  introduced  into  the  mixture,  on  the 
basis  of  a  specialized  checking  procedure,  such  that  it 
is  the  condensed  species  not  included  in  the  mixture 
that  should  lower  L(n,  A)  the  most  if  added  to  the  mix¬ 
ture.  Then  T(n,  A)  is  minimized  again  with  the  addi¬ 
tional  condensed  species  (initial  mole  number  of  zero), 
to  perform  the  second  iterative  loop.  If  any  condensed 
species  included  in  the  mixture  has  a  negative  mole 
number,  it  is  removed  from  the  mixture.  If  another 
condensed  species  not  included  in  the  mixture  should 
be  included,  this  is  done,  and  another  iterative  loop  is 
performed.  This  looping  ends  and  the  final  equilibrium 
composition  n  is  obtained  when  condensed  species  are 
neither  removed  nor  added  to  the  mixture. 

The  iterative  procedure  in  BNR  methods  is  an 
inefficient  multiple-loop  structure  that  cannot  be  re¬ 
duced  to  a  single  loop,  because  the  specialized  checking 
procedure  for  condensed  species  is  based  on  using  equi¬ 
librium  mixture  conditions  that  are  established  only  at 
the  end  of  an  iterative  loop.  In  the  NASA-CEA  and 
Bagheera  codes  the  variables  In(ni)  are  used  instead  of 
rii  for  the  gaseous  species  only,  which  has  the  impor¬ 
tant  advantage  of  making  all  gaseous  mole  numbers 
nonnegative  during  the  iterative  procedure. 

The  CERV  method  and  code  determine  the  equi¬ 
librium  mixture  composition  for  a  T-P  problem  by 
minimizing  the  Gibbs  energy  G(^)  expressed  in  terms 
of  reaction  variables  or  stoichiometric-vector  multipli¬ 
ers  numbering  =  Ng  —  Ng  instead  of  Ag-l-iVe  com¬ 
position  variables  and  Lagrange  constraint  multipliers. 
The  reaction  variables  are  related  to  the  composition 
variables  by  the  linear  transformation  n  =  n*  N^, 
where  n*  is  a  known  set  of  mole  numbers  (e.g.,  initial 
guesses)  and  Af  is  a  special  Ng  x  matrix  of  stoichio¬ 
metric  coefficients.  The  mass  conservation  constraints 
on  the  atomic  elements  {Fn  =  a)  are  integrated  di¬ 
rectly  into  the  formulation  by  means  of  matrix  N,  and 
each  column  of  TV  is  a  set  of  stoichiometric  coefficients 
associated  with  a  particular  chemical  reaction  occur¬ 
ring  among  some  gaseous  and  condensed  species. 

The  minimum  of  G(^)  and  corresponding  equilib¬ 
rium  composition  n  are  obtained  iteratively  by  using  a 
standard  second-order  variation  scheme.  An  appropri¬ 
ate  set  of  gaseous,  liquid  and  solid  species  is  selected 
as  reaction  products,  and  their  mole  numbers  are  esti¬ 
mated  by  using  linearization  and  the  well-known  sim¬ 
plex  method  that  ensures  the  mass  conservation  con¬ 


straints  are  satisfied.  A  quadratic  approximation  Q(4) 
to  G(4)  is  constructed  in  terms  of  by  matching  the 
two  functions  in  terms  of  n  and  first  derivatives  of  G(4) 
with  respect  to  Finding  the  minimum  of  Q{$,)  leads 
to  a  set  of  TVr  linear  equations  of  the  form  HA^  =  c, 
in  which  JT  is  a  TV^  x  iVr  Hessian  matrix.  In  the  CERV 
method  the  Hessian  matrix  is  made  diagonally  dom¬ 
inant  because  of  the  special  form  of  matrix  TV.  The 
off-diagonal  elements  are  neglected,  the  approximate 
inverse  is  obtained  easily  and  efficiently  by  in¬ 
verting  TVr  scalar  and  nonzero  diagonal  numbers,  the 
solutions  =  H~^c  and  An  =  TVA^  are  calcu¬ 
lated,  and  then  ^  -I-  A^  and  n  =  n*  An  are 

obtained  for  the  minimum  of  Q{0-  The  new  values  ^ 
and  n  are  closer  than  the  old  ones  and  n*  to  those 
that  define  the  minimum  of  G(^).  Repetition  of  the 
quadratic  procedure  produces  a  step-by-step  descent 
path  to  the  minimum  of  G(4)  and  a  converged  solu¬ 
tion  for  the  equilibrium  composition  n.  See  the  book 
by  Smith  and  Missen  (1991)  and  the  report  by  Wong, 
Gottlieb  and  Lussier  (2003)  for  more  details  on  the 
minimization  procedure  and  CERV  method. 

Owing  to  the  reaction-variable  formulation  of  the 
CERV  method,  the  associated  computer  code  does  not 
suffer  convergence  failures.  The  method  also  features 


Table  2.  Comparison  of  BNR  and  CERV 
methods  for  T-P  problems. 


feature  \  method 

BNR 

CERV 

solution  method 

minimize  L{n,  A) 
by  Newton  Raph- 
son  scheme 

minimize  G(^) 
by  second  var- 
-iation  scheme 

variables 

n  or  ln(n),  A 

n 

equations 

As-b  Ae 

Ae  -  Ae 

iteration  type 

multiple  loop 

single  loop 

mass  conservation 
of  atomic  elements 
{Fn  =  a) 

side 

constraint 
using  A 

embedded 
constraint  in 
n  =  n*  -b 

mass  conservation 
satisfied  every  step 

no 

yes 

equilibrium  equation 

-  F^\  =  0 

TV^/x  =  0 

equilibrium  equation 
satisfied  every  step 

yes  —  Brinkley 
no  —  NASA 
no  —  RAND 

no 

nonnegativity 

constraint 

handled 
very  easily 

handled 

easily 

chemical  reaction 
equations  included 

no 

yes 

add  and  remove 
condensed  species 

once  per  loop, 
reconstruct  F 

once  per  step, 
reconstruct  TV 

theoretically  singular 
matrix  that  causes 
convergence  failure 

must  avoid 
at  outset 

NA 

numerically  singular 
matrix  that  causes 
convergence  failure 

possible 

NA 

convergence  enforcer 

yes 

yes 

computational  speed 

good 

faster 
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an  efficient  single-loop  iterative  structure,  because  con¬ 
densed  species  can  be  removed  from  or  added  to  the 
reaction  products  at  each  step  along  the  single  de¬ 
scent  path,  because  the  mass  constraints  are  satisfied 
at  each  step  of  the  single  loop  and  because  the  special¬ 
ized  checking  procedure  is  based  on  species  chemical 
reactions  that  are  embedded  in  matrix  TV  (rather  than 
on  equilibrium  conditions  as  in  BNR  methods). 

The  essential  features  of  BNR  and  CERV  methods 
for  solving  T-P  problems  are  given  in  Table  2.  These 
results  highlight  method  similarities  and  differences, 
and  some  major  advantages  of  the  CERV  method  and 
code.  The  robustness  of  the  CERV  code  in  not  suf¬ 
fering  convergence  failure  is  especially  important  in 
military  applications  involving  primers,  fuses  and  ig¬ 
niters.  The  single-loop  iterative  structure  combined 
with  the  simple  and  quick  Hessian  matrix  inversion 
result  in  the  most  efficient  computations  for  equilib¬ 
rium  reaction  products  modelled  as  nonideal  solutions 
of  imperfect  gases  and  compressible  condensed  species. 

4.  COMPUTATIONAL  RESULTS 
4.1  E-V  Problem  for  PETN 

The  problem  of  the  detonation  of  pentaerythri- 
tol  tetranitrate  (PETN)  in  a  closed  vessel  with  atmo¬ 
spheric  air  (298.15  K,  101.325  kPa)  surrounding  the 
charge  is  solved  by  the  CERV  code.  The  final  equi¬ 
librium  temperature,  pressure  and  mixture  composi¬ 
tion  are  determined  iteratively  by  repeatedly  solving 
a  T-P  problem  and  systematically  adjusting  the  mix¬ 
ture  temperature  while  updating  the  mixture  pressure 
based  on  the  known  closed-vessel  volume  at  each  itera¬ 
tive  step  until  the  final  energy  of  the  reaction  products 
equals  the  known  initial  energy  of  the  three  reactants 
C5H8N4O12,  N2  and  O2.  Other  ingredients  of  air  such 
as  Ar  and  CO2  are  small  and  ignored,  and  the  ignition 
source  for  PETN  is  neglected. 

The  computations  are  done  for  closed-vessel  load¬ 
ing  densities  S  varying  from  0.01  to  0.20  gm/cm^, 
where  6  denotes  the  mass  of  PETN  divided  by  the 
vessel  volume.  Results  are  shown  in  Figs.  1  and  2  for 
the  cases  of  a  nonideal  mixture  model  and  four  dif¬ 
ferent  equations  of  state.  When  the  nonideal  mixture 
model  is  combined  with  the  simplest  EOS  given  by 
Pv  =  TZT  for  perfect  gases  (PC),  the  model  for  the 
nonideal  mixture  reduces  to  that  for  an  ideal  mixture 
of  perfect  gases,  and  this  serves  as  a  good  reference. 

The  equilibrium  mixture  pressure  P  for  the  case  of 
a  perfect  gas  increases  almost  linearly  with  increasing 
vessel  loading  density,  as  shown  in  Fig.  1,  illustrating 
that  T/v  also  increases  linearly  with  loading  density. 
The  pressure  variations  for  the  cases  of  the  other  EOS 
all  rise  more  rapidly;  this  steeper  rise  stems  primarily 
from  the  effects  of  co-volume  (stronger  intermolecular 
repulsion  forces  for  larger  gas  densities).  The  equi¬ 
librium  mixture  temperature  T  shown  in  Fig.  1  fea¬ 
tures  an  initial  rapid  increase  for  loading  densities  be- 


Figure  1:  Equilibrium  temperature  and  pressure 
for  the  detonation  of  PETN  in  a  closed  vessel. 


Figure  2:  Equilibrium  mole  fractions  for 
the  detonation  of  PETN  in  a  closed  vessel. 


low  0.05  gm/cm^,  but  this  initial  rapid  rise  levels  off 
quickly  and  appears  to  approach  a  maximum  temper¬ 
ature  asymptotically.  The  temperature  should  trend 
toward  a  constant  adiabatic  flame  temperature  as  the 
loading  density  increases,  or  as  this  temperature  be¬ 
comes  less  sensitive  to  an  increasing  loading  density. 

The  most  abundant  reaction  products  for  this 
closed- vessel  problem  are  the  gaseous  species  H2O, 
CO2,  CO,  N2,  OH,  H2,  O2,  NO,  H,  and  O,  and  their 
mole  fractions  are  shown  in  Fig. 2,  although  many  other 
gaseous  species  such  as  CH4,  NO2,  N2O  and  N  are  in¬ 
cluded  in  the  computations.  For  the  combination  of 
atomic  elements  C,  H,  N  and  O  of  the  reactants,  the 
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possible  condensed  species  are  C(s-gr),  H20(liquid) 
and  H2O (solid).  The  abundance  of  graphite  was  negli¬ 
gible,  and  liquid  water  and  ice  do  not  exist  at  the  high 
equilibrium-mixture  temperatures  for  this  problem. 

The  species  mole  fractions  shown  in  Fig.  2  do  not 
appear  to  be  influenced  greatly  by  the  four  different 
EOS,  but  nontrivial  changes  are  obscurred  by  the  log¬ 
arithmic  scale.  The  maximum  differences  are  2.3,  2.7, 
2.0,  0.5  and  31%  among  the  H2O,  CO2,  CO,  N2  and 
O2  data,  respectively,  at  5  =  0.20gm/cm^. 

4.2  T-P  Problem  with  Solid  or  Liquid  Species 

A  T-P  problem  is  formulated  from  the  reactants 
that  are  composed  of  the  four  atomic  elements  C,  K, 
O  and  S.  The  reactants  are  taken  as  the  ten  gaseous 
species  CO2,  CO,  SO2,  K,  SO,  O2,  S,  C,  K2SO4  and 
K2,  each  in  the  amount  of  one  mole,  and  the  four 
condensed  species  as  K2S04(s),  K2S04(lq),  K2C03(s) 
and  K2C03(lq),  each  in  the  amount  of  one  tenth  of 
a  mole.  Equilibrium-mixture  compositions  are  com¬ 
puted  by  the  CERV  code  for  a  specified  mixture  pres¬ 
sure  P  =  100  kPa  (1  bar)  and  two  specified  mixture 
temperatures  T  =  700  and  1500  K.  The  two  sets  of 
computed  results  are  shown  in  the  bar  chart  in  Fig.  3, 
in  which  the  gaseous  species  are  grouped  at  the  top  and 
the  condensed  species  are  assembled  at  the  bottom. 
Note  that  virtually  identical  results  are  obtained  for 
computations  based  on  an  ideal  or  a  nonideal  mixture 
of  perfect  or  imperfect  gases  because  the  equilibrium- 
mixture  pressure  is  relatively  low  for  this  problem. 

The  gaseous  species  CO,  S,  SO,  CO2  and  SO2 
occur  in  abundance  larger  than  10“®  moles  at  both 
temperatures  of  700  and  1500  K.  However,  at  700  K 
the  amounts  of  the  gaseous  species  K,  K2SO4  and 
K2  are  smaller  than  10“®  moles,  and  the  two  sub¬ 
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Figure  3:  Mole  numbers  of  gaseous  and  condensed 
species  for  a  mixture  at  P  =  100  kPa  and  T  =  700 
and  1500  K. 


stances  potassium  sulphate  K2SO4  and  potassium  car¬ 
bonate  K2CO3  occur  as  solid  species  in  large  amounts, 
whereas  their  liquid  counterparts  are  nonexistent.  At 
1500  K  the  gaseous  species  K,  K2SO4  and  K2  occur 
in  more  significant  amounts  and  the  two  substances 
K2SO4  and  K2CO3  now  exist  as  liquid  species  only. 
This  illustrates  that  large  changes  in  the  equilibrium 
mixture  temperature  can  result  in  large  changes  in  the 
amounts  of  some  gaseous  species,  and  some  substances 
can  change  in  form  from  solid  to  liquid  species. 

For  this  particular  T-P  problem  the  species  in  the 
reaction  products  are  arbitrarily  taken  as  identical  to 
those  in  the  reactants.  The  mole  numbers  of  all  re¬ 
actant  species  are  specified  at  the  outset.  The  initial 
set  of  mole  numbers  of  the  reaction-product  species 
can  be  estimated  as  mentioned  previously  by  using 
the  linearization  and  simplex  method,  or  they  can  be 
taken  simply  as  equal  to  those  of  the  reactants  because 
the  reactant  and  reaction-product  species  are  identical. 
This  T-P  problem  is  solved  both  ways  by  the  CERV 
code,  and  the  final  results  for  the  two  temperatures 
shown  in  Fig.  3  are  the  same,  as  might  be  expected. 

This  exercise  helps  demonstrate  that  the  CERV 
method  and  code  can  handle  multiple  species  with  the 
same  chemical  formula  (i.e.,  liquid  and  solid  forms  of 
substances  K2SO4  and  K2CO3)  as  part  of  the  reaction 
products,  and  only  the  naturally  occurring  form  of  a 
substance  will  exist  in  the  final  equilibrium  mixture. 
This  type  of  input  cannot  be  used  by  codes  based  on 
the  BNR  method,  because  matrix  A  would  be  theo¬ 
retically  singular  and  produce  a  convergence  failure. 

4.3  E-V  Problem  with  Solid-Liquid  Transition 

The  black  powder  for  this  closed-vessel  problem 
is  taken  as  composed  by  weight  of  45%  solid  potas¬ 
sium  nitrate  (KNO3),  23%  solid  Susquehanna  charcoal 
(C1321H573S6N4O100),  17%  solid  sulphur  (S),  0.5%  liq¬ 
uid  water  (H2O)  and  14.5%  solid  ash  or  potassium  car¬ 
bonate  (K2CO3)  at  298. 15K  and  101.325 kPa.  The  re¬ 
actants  with  the  six  atomic  elements  C,  H,  K,  N,  O  and 
S  involve  156  gaseous  species,  29  condensed  species  and 
12  pairs  with  possible  phase  transitions  between  con¬ 
densed  species  of  the  same  substance.  For  demonstra¬ 
tion  we  select  the  melting  of  solid  potassium  sulphide 
K2S(s-3)  to  liquid  K2S(lq)  as  the  phase  transition  pair, 
with  a  melting  temperature  Tmeit  =  1222.24  K. 

The  equilibrium-mixture  temperature  T,  pressure 
P  and  composition  n  are  computed  by  the  CERV 
code  for  the  closed-vessel  loading  density  5,  mass  of 
the  black  powder  divided  by  the  vessel  volume,  vary¬ 
ing  from  0.096  to  0.105  gm/cm^.  Some  computed  re¬ 
sults  are  presented  in  Fig.  4.  For  loading  densities 
from  0.096  to  0.09726  gm/cm^  the  mixture  pressure 
and  temperature  increase,  and  only  the  solid  phase 
K2S(s-3)  of  potassium  sulphide  occurs  in  the  mix¬ 
ture.  For  0.09726  <  6  <  0.10377  gm/cm^  the  mixture 
temperature  remains  fixed  at  the  melting  temperature 
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Figure  4:  Combustion  of  black  powder  in  a  closed 
vessel  with  the  melting  of  potassium  sulphide. 


1222.24  K,  the  mixture  pressure  continues  to  increase, 
the  amount  of  the  solid  phase  K2S(s-3)  diminishes  to¬ 
ward  zero,  and  the  amount  of  the  liquid  phase  K2S(lq) 
increases  from  zero  as  melting  continues.  The  solid 
phase  disappears  at  the  end  of  this  range  and  only  the 
liquid  phase  exists  afterwards.  For  loading  densities 
larger  than  0.10377  gm/cm^  the  mixture  temperature 
increases  once  again,  along  with  the  mixture  pressure. 

4.4  T-P  Problem  with  Condensation 

Consider  a  mixture  of  inert  nitrogen  N2,  benzene 
CeHg  and  toluene  CrHg  at  a  specified  mixture  tem¬ 
perature  and  range  of  mixture  pressures  such  that  the 
benzene  and  toluene  can  occur  in  gaseous  and  liquid 
forms.  The  reactants  are  one-half  mole  of  gaseous 
nitrogen  and  one  mole  each  of  gaseous  benzene  and 
toluene.  The  mixture  composition  n  is  calculated 
by  the  CERV  code  for  the  reaction  products  N2(g), 
C6H6(g),  C7H8(g),  C6H6(lq)  and  C7H8(lq),  for  a  speci¬ 
fied  mixture  temperature  T  =  298. 15K  and  pressure  P 
ranging  from  5  to  30kPa.  Computed  results  are  shown 
in  Fig.  5,  along  with  the  partial  pressures  PceHaig)  E^nd 
PCjHg{g)  of  gaseous  benzene  and  toluene.  The  vapour 
pressures  12.7  and  3.78  kPa  of  benzene  and  toluene, 
respectively,  are  also  indicated  in  the  figure. 

Only  gaseous  nitrogen,  benzene  and  toluene  occur 
in  the  equilibrium  mixture  at  mixture  pressures  be¬ 
low  9.50  kPa.  At  the  mixture  pressure  of  9.50  kPa  the 
gaseous  toluene  begins  to  condense  into  liquid  toluene, 
and  this  condensation  continues  with  increasing  mix¬ 
ture  pressure,  resulting  in  a  decrease  in  the  abundance 
of  gaseous  C7H8  and  a  corresponding  increase  in  the 
amount  of  liquid  C7H8.  At  the  mixture  pressure  of 
22.85  kPa  the  condensation  of  gaseous  benzene  into 
liquid  benzene  begins,  and  this  condensation  contin¬ 
ues  as  the  mixture  pressure  increases. 

When  the  chemical  potential  of  the 

gaseous  species,  calculated  on  the  basis  of  the  gaseous 


mixture  pressure  P  (kPa) 


Figure  5:  Condensation  of  benzene  and  toluene. 

mole  fraction  Xi  and  mixture  temperature  T  and  pres¬ 
sure  P,  is  less  than  the  chemical  potential  p^{p'j,  T,  P) 
of  the  related  liquid  species,  calculated  on  the  basis  of 
its  vapour  pressure  pj  and  mixture  temperature  T  and 
pressure  P,  then  only  the  gaseous  species  exist  in  the 
mixture.  When  these  two  chemical  potentials  just  be¬ 
come  equal,  as  the  mixture  pressure  increases,  conden¬ 
sation  just  begins.  For  higher  mixture  pressures  the 
condensation  continues,  and  the  relative  amounts  of 
these  gaseous  and  liquid  species  are  determined  partly 
on  the  basis  of  p,f{xi,P)  =  p*'^(pJ,P).  For  an  ideal 
mixture  of  perfect  gases  and  a  simplified  liquid  model, 
these  general  relationships  reduce  to  using  the  partial 
pressure  of  the  gaseous  species  pi  =  XiP  instead  of 
p.f{xi,  P)  and  the  vapour  pressure  of  the  liquid  pj  in¬ 
stead  of  P^{p"ji  P)-  For  this  problem  with  low  mixture 
pressures  and  densities,  these  simplified  results  are  suf¬ 
ficient  to  explain  the  results  shown  in  Fig.  5. 

4.5  E-V  Problem  with  Compressible  Species 

The  black  powder  in  the  closed  vessel  is  composed 
by  weight  of  69%  solid  potassium  nitrate  KNO3,  25% 
solid  graphite  C(s-gr)  and  6%  solid  sulphur  S.  These 
reactants  are  initially  surrounded  by  air  (79%  N2  and 
21%  O2)  at  298.15  K  and  101.325  kPa.  The  reac¬ 
tion  products  composed  of  the  six  atomic  elements 
C,  H,  K,  N,  O  and  S  once  again  involve  156  gaseous 
species,  29  condensed  species  and  12  pairs  with  possi¬ 
ble  phase  transitions  between  condensed  species  of  the 
same  substance.  The  most  abundant  gaseous  species 
in  this  problem  are  CO,  N2,  CO2  KCN,  K,  COS, 
K2C2N2,  K2CO3,  CS2,  K2,  S2,  CS,  SO2,  SO  and 
S2O.  These  gaseous  reaction  products  are  modelled 
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as  a  nonideal  mixture  of  imperfect  gases  by  using  a 
three-term  truncated  virial  equation  of  state.  The  con¬ 
densed  species  that  occur  in  the  equilibrium  mixture 
are  solid  graphite  C(s-gr),  liquid  potassium  carbonate 
K2C03(lq)  and  liquid  potassium  sulphide  K2S(lq). 

Compressibility  effects  of  condensed  species  are  il¬ 
lustrated  by  CERV  computations  with  three  compress¬ 
ibility  models.  The  simplest  no-volume  model  (used  in 
the  NASA-CEA  code)  includes  the  mass  and  energy  of 
the  condensed  species  but  neglects  their  volume  (i.e., 
material  density  is  infinite).  The  constant-volume  or 
incompressibility  model  (used  in  the  Bagheera  code) 
includes  the  mass  and  energy  of  the  condensed  species 
with  a  constant  volume  and  density  evaluated  for  ei¬ 
ther  the  solid  or  liquid  that  exists  at  the  initial  tem¬ 
perature  and  pressure  of  the  reactants  in  the  closed 
vessel.  The  variable-volume  or  compressibility  model 
(used  in  the  Blake  code)  incorporates  an  approximate 
pressure-volume-temperature  equation  of  state  to  in¬ 
clude  the  counter-acting  effects  of  thermal  expansion 
and  pressurized  contraction  on  the  condensed  species 
at  high  mixture  temperatures  and  pressures.  Equa¬ 
tions  of  state  for  compressible  condensed  species  were 
adopted  from  the  Blake-code  library. 

Computational  results  for  the  equilibrium-mixture 
temperature  and  pressure  are  presented  in  Fig.  6.  The 
lowest  mixture  temperatures  and  pressures  occur  for 
condensed  species  with  no  volume,  because  this  leaves 
more  volume  for  the  gaseous  species  to  occupy  at  a 
lower  mixture  pressure.  For  const  ant- volume  or  in¬ 
compressible  condensed  species  the  mixture  tempera¬ 
ture  and  pressure  increase  by  0.2%  and  5.9%,  respec¬ 
tively,  at  (5  =  0.20  gm/cm^.  The  condensed-species 
volume  reduces  the  volume  available  to  the  gaseous 
species,  which  has  only  a  small  effect  on  the  mix¬ 
ture  temperature  but  an  impact  on  the  mixture  pres- 


Figure  6:  Volume  effects  of  condensed  species  on 
equilibrium-mixture  temperatures  and  pressures. 


Figure  7:  Volume  effects  of  condensed  species  on 
equilibrium-mixture  species  abundances. 


sure.  For  variable- volume  or  compressible  condensed 
species  the  mixture  temperature  and  pressure  are  the 
highest,  being  0.3%  and  7.4%  larger  than  for  the  no¬ 
volume  model.  The  increase  in  mixture  pressure  from 
changing  from  constant-  to  variable-volume  condensed 
species  stems  mainly  from  the  large  thermal  expan¬ 
sion  of  species  material  at  high  mixture  temperatures, 
which  is  only  partly  off-set  by  the  species-material 
compression  at  high  mixture  pressures,  yielding  overall 
a  smaller  volume  available  to  the  gaseous  species  that 
gives  a  higher  mixture  pressure.  The  mixture  temper¬ 
ature  is  insensitive  to  the  three  compressibility  models 
because  the  mass  and  energy  of  the  condensed  species 
are  included  in  all  three  models. 

The  effects  of  compressibility  on  the  gaseous  and 
condensed-species  abundances  are  presented  in  Fig.  7. 
The  maximum  changes  in  abundances  among  the 
three  compressibility  models  at  a  loading  density  of 
0.20gm/cm^  are  1.8%  for  CO2,  0.7%  for  K2S(lq),  0.5% 
for  K2C03(lq),  0.4%  for  CO  and  C(s-gr),  7.6%  for 
COS  and  0%  for  N2.  These  effects  on  the  most  abun¬ 
dant  gaseous  and  condensed  species  are  therefore  quite 
small.  Note  that  carbon  monoxide  (CO)  is  the  most 
abundant  species  by  a  factor  of  about  two  to  six,  and 
its  abundance  only  is  given  by  the  larger  scale. 

5.  CONCLUDING  REMARKS 

Major  computational  codes  called  NASA-CEA, 
Blake,  Bagheera  and  CERV  and  their  solution  meth¬ 
ods  were  described  to  highlight  their  capabilities  in 
solving  civilian  and  military  problems  associated  with 
engineering  and  energetic-materials  applications,  and 
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also  to  illustrate  that  modern  codes  need  to  solve 
closed-vessel  problems  without  convergence  failure 
with  hundreds  of  gaseous  and  condensed  species  using 
models  of  nonideal  mixtures  of  imperfect  gases,  com¬ 
pressible  condensed  species  and  liquid-solid  and  solid- 
solid  phase  transitions.  The  CERV  method  and  code, 
based  on  reaction  instead  of  composition  variables, 
has  inherent  advantages  mentioned  earlier  in  achiev¬ 
ing  these  goals.  Although  good  methods  and  codes  are 
available,  computations  are  nonetheless  limited  by  the 
scarcity  of  data  for  gaseous  species  (EOS  for  imperfect 
gases,  mixing  rules)  and  condensed  species  (EOS  for 
thermal  expansion  and  pressurized  contraction). 

Some  closed- vessel  or  E-V  problems  were  intro¬ 
duced  and  solved  so  that  results  could  be  presented 
to  illustrate  the  effects  of  nonideal  mixtures  of  perfect 
and  imperfect  gases  and  the  effects  of  compressibility 
of  condensed  species  on  the  equilibrium-mixture  tem¬ 
peratures,  pressures  and  compositions.  The  effects  of 
compressible  in  contrast  to  incompressible  condensed 
species  on  enhancing  the  mixture  pressure  from  the 
predominant  effect  of  thermal  expansion  over  high- 
pressure  material  compression  is  interesting  and  re¬ 
ported  herein  for  the  first  time.  However,  more  data 
is  required  to  address  these  issues  more  completely. 

Two  T-P  problems  were  introduced  to  illustrate 
that  the  CERV  method  and  code  can  simultaneously 
handle  gaseous,  liquid  and  solid  species  of  the  same 
substance  with  ease,  whether  a  particular  problem  has 
a  final  equilibrium  mixture  that  features  only  the  solid 
species,  only  the  liquid  species,  both  gaseous  and  liquid 
species  undergoing  condensation,  both  liquid  and  solid 
species  in  a  phase  transition  (melting),  or  two  solid 
species  in  a  phase  transition  (not  demonstrated) . 

The  capability  of  a  computer  code  to  readily  solv¬ 
ing  T-P  problems  with  multiple  gaseous,  liquid  and 
solid  species  formed  from  the  same  substance  is  impor¬ 
tant  to  the  solution  of  E-V  or  closed-vessel  problems. 
The  solution  of  an  E-V  problem  is  determined  itera¬ 
tively  by  repeatedly  solving  a  T-P  problem  and  sys¬ 
tematically  adjusting  the  mixture  temperature  while 
updating  the  mixture  pressure  based  on  the  known 
closed-vessel  volume  at  each  iterative  step  until  the 
final  energy  of  the  reaction  products  is  equal  to  the 
initial  energy  of  the  reactants.  Each  new  T-P  prob¬ 
lem  should  use  the  mole  numbers  of  the  previous  T-P 
problem  as  initial  guesses  for  computational  efficiency 
during  the  iterative  procedure.  Because  a  previous  so¬ 
lution  may  contain  multiple  gaseous,  liquid  and  solid 
species  formed  from  the  same  substance  during  con¬ 
densation  or  melting,  they  need  to  be  handled  without 
difficulty  or  convergence  failure  in  proceeding  to  the 
next  iterative  step.  This  is  achieved  easily  by  means 
of  the  CERV  methodology  (with  reaction  variables) 
and  with  a  well- written  CERV  code. 

The  BNR  methodology  (with  composition  vari¬ 
ables)  in  contrast  has  inherent  difficulties  in  using  mole 


numbers  directly  from  an  earlier  iterative  step  when 
multiple  liquid  and  solid  species  of  the  same  substance 
are  present,  because  matrix  A  can  become  singular 
theoretically  or  numerically  with  a  convergence  failure 
on  inversion,  unless  some  species  are  rejected  from  the 
mixture  to  make  A  invertible  before  proceeding. 
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